Towards a first-principles theory of surface thermodynamics and kinetics 
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Understanding of the complex behavior of particles at surfaces requires detailed knowledge of both 
macroscopic and microscopic processes that take place; also certain processes depend critically on 
temperature and gas pressure. To link these processes we combine state-of-the-art microscopic, 
and macroscopic phenomenological, theories. We apply our theory to the O/Ru(0001) system 
and calculate thermal desorption spectra, heat of adsorption, and the surface phase diagram. The 
agreement with experiment provides validity for our approach which thus identifies the way for a 
predictive simulation of surface thermodynamics and kinetics. 

PACSnumbers: 68.45. Da, 82.65. My, 82.65.Dp 



The study of the physical and chemical processes that 
take place at gas-surface interfaces have long been an area 
of intense research. This interest is both fundamental 
as well as driven by the possible discovery of important 
technological applications, e.g.p.iii the field of heteroge- 
neous catalysis, corrosion, etcBEl. With respect to the 
field of the theory of adsorption of gases on solid sur- 
faces, advancement in recent years has developed in two 
distinct, albeit complimentary directions: (i) electronic 
structure calculations, at best done by density-functional 
theory (DPT), to determine the geometries, energetics, 
and vibrational properties of adsorbate covered surfaces, 
and (ii) phenomenologicaL-podels, both for the thermo- 
dynamics and the kinetica3 of the adsorbate. If one can 
assume that the geometry of the solid surface does not 
change dramatically and that adsorption occurs at well 
defined sites, one frequently employs a lattice gas model. 
A number of parameters enter such a model such as the 
binding energies and vibrational frequencies of a single 
adparticle in the various adsorption sites, and their mu- 
tual lateral interactions with adparticles in close-by sites. 
Traditionally, these parameters are adjusted in the the- 
ory to fit a variety of experimental data such as phase 
diagrams, heats of adsorption, infrared spectra and ther- 
mal desorption data, etc. Such an approach, while use- 
ful, is clearly not necessarily predictive in nature, nor the 
parameters unique, and may not capture the physics of 
the microscopic processes that are behind the "best-fit" 
adjusted "effective" parameters. 

In the present letter, with the aim to improve upon 
this approach, we combine state-of-the-art procedures of 
(i) microscopic theories, i.e. DPT electronic structure 
calculations and (ii) macroscopic phenomenological ap- 
proaches, i.e. lattice gas and rate equations, and Monte 
Carlo schemes. On doing this, we present a consistent 
first- principles based approach for calculation of the ther- 
modynamic and kinetic properties of an adsorbate, such 
as heats of adsorption, temperature programmed des- 
orption (TPD) spectra, and the surface phase diagram. 
We have chosen the systjera. of oxygen at Rn4p001) for 
which detaiJjg^ structural B, thermodynamidlj, and ki- 



We like to mention that use of results of ab initio elec- 
tronic structure calculations as input to statistical me- 
chanical computer codes have also been dpzjeloped for the 
description of bulk alloy pbase diagramstsEj and island 
growth on crystal surfacesO. To our knowledge, how- 
ever, there has been no application of such an approach 
for the processes that we are interestedjn here. 

The electronic structure calculationsllj are performed 
using DPT and the generalized gradient approppiation 
(GGA) for the exchange-correlation functionalllll (here- 
aftee, denoted as DPT-GGA). We use the pseudopoten- 
tiallaH plane wave method and the supercell approach 
to model the surface. The position of the O atoms and 
the top two Ru layers are fully relaxed. The DPT-GGA 
Ru pseudopotential yields a bulk hcp-fcc energy differ- 
ence of « — (iX)72 eV in good agreement with all-clcfjtijon 
calculations .Ej Por further details we refer to Refs.HEj 

To set up a lattice gas model for O on Ru (0001) we 
require a hamiltonian, which we express as: 
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netic dataUO exist. We will show that with the present 
approach, a realistic description of these physical prop- 
erties is indeed feasible. 
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The different unit cells of the substrate surface are la- 
beled by an index i with i -\- a, etc. labeling neighboring 
cells, and we introduce occupation numbers = or 1 
depending on whether a site in cell i is empty or occu- 
pied. Overcounting of cells is excluded in the summa- 
tions. Eq. 1 includes consideration of hep and fee sites. 
The indicies a', 6', and c' indicate that the first, second, 
and third neighbor distances between atoms in hep and 
fee sites are different to when they occupy the same type 
of sites. Here E^'^p — \Vo\ — ksTinq^ is the binding en- 
ergy of an isolated particle in an hep site. |Vb| is the 
depth of the adsorption potential with reference to the 
energy of a gas phase molecule which adsorbs dissocia- 
tively. qs is the partition function of the atom accounting 
for its vibration perpendicular, and its frustrated trans- 



lation parallel, to the surface. V^^^, V2n^ and V^^^ are 
the first, second and third neighbor interaction energies 
between two adsorbed O atoms in hep sites. The anal- 
ogous terms labeled "fee" represent the same quantities 
but for adsorption in fee sites. Terms labeled ]/hcp-fcc 
represent the interactions between atoms in hep and fee 
sites. Trio (and higher cluster) interaction energies, Vtrio, 
account for additional modifications because the interac- 
tion between two adsorbed O atoms is changed, when a 
third adatom is close by. In fact, depending on the angles 
and distances, there are three different trio interactions 
taken into account. 

We determine the lateral interaction energies required 
in Eq. 1 from DFT-GGA calculations of ordered struc- 
tures of O on Ru(OOOl) (see Fig. 1). The adsorption 
energy of a single oxygen atom on Ru (0001) is obtained 
using a (3 X 3) structure with coverage 9 = 1/9. With 
this large 0-0 separation (8.26 A), which corresponds to 
the fifth nearest neighbor distance of alike sites, lateral 
interactions are negligible. The adsorption energy, with 
respect to I/2O2, is then expressed as, 

j^e=l/9 _ pO/Ru _ pRu 
^" — -^total = 



^total 
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Here £^^(2", £;£tai> and E^J^^^ are the total energies of 
the (3x3)-O/Ru(0001) system, the clean Ru surface, and 
a free O2 molecule, respectively. Expressions analogous 
to Eq. 2 have been used to derive the adsorption ener- 
gies for the other coverages (see Tab. I). We expand the 
adsorption energies in terms of two- and three-body in- 
teractions. The interaction parameters are derived from 
the equations: 
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= Vo + Vm + V2n + 3^3„ + Vlt 



=Vo + IVm + 3F2n + I Fan + 31/bt 
=3/4 = Vo+ 2Vin + 2V2„ + 3F3n + 2Vlt + 2Vbt + iVtt 

Et^° = Vo + HVm + V2n + V3n) + Wit + e^bt + 2Vtt 

where Vn, Vbt and Vtt are linear, bent and triangular 
trios, respectively (indicated in Fig. 1). We point out 
that only in the (3 x 3)-20 structure at 6* = 2/9 can the 
atoms at nearest neighbor sites move significantly from 
the locally 3-fold symmetric adsorption sites to reduce 
the repulsion, Vin- However, with strong nearest neigh- 
bor repulsion, isolated nearest neighbor pairs are highly 
improbable at any coverage. 



We have therefore calculated Ea~^^^ for atoms at lo- 
cally 3-fold symmetric sites because if we were to use the 
relaxed structure we would need to include higher many- 
body interactions (e.g. longer-ranged trios, quartos and 
quintos etc.) to account for the movement of the atoms 
back to the ideal 3-fold sites which occurs for higher cov- 
erages. 

We thus have for O atoms in hep sites (and analo- 
gously for O atoms in fee sites) seven equations for six 
unknowns. Using the first six of these equations we ob- 
tain the interaction energies listed in Tab. II. The accu- 
racy of these interaction values is gauged by calculating 
the adsorption energy of the monolayer (seventh equa- 
tion); we find a discrepancy of 0.034 eV smaller than 
that of the DFT-GGA value. Along with our derived 
values in Tab. II, we give in brackets the interaction en- 
ergies determined by Piercy et al&3 for their best-fit to 
the experimental O/Ru(0001) phase diagram. We also 
include our derived interaction parameters for structures 
involving O atoms in hep and fee sites (last three dia- 
grams of Fig. 1) obtained by writing down appropriate 
equations in an analogous manner to those listed above. 

To complete the specification of our hamiltonian, we 
need the vibrational frequencies of an O atom relative 
to the Ru surface. These frequencies can also be calcu- 
lated using density-functional theory. For example, we 
obtain for the vibration of oxygen normal to the surface 
(at the F-point) i^z = 509 cm~^. The experimental value 
is 535 cm-i (for O in the (2 x 2) phase)M 

We now proceed to calculate the temperature pro- 
grammed thermal desorption spectra. Writing the ki- 
netic equation for adsorption and desorption as dO/dt = 
-Rad— Rdcs, we obtain for an atomic adsorbate in contact 
with a gas of diatomic homonuclear molecules, the rate of 
adsorption i?ad — 2Sdis{9,T)PmasXm/h. Here Pm is the 
molecular pressure above the surface, as is the area of one 
unit cell of the substrate surface, Am = ft,/(27rmA:sT)^/^ 
is the thermal wavelength of a molecule of mass m, and 
Sdis{0,T) is the dissociative sticking coefficient. For the 
rate of desorption we haveB, 
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Here Z^r is the partition function accounting for the in- 
ternal vibrations and the rotations of O2 in the gas phase. 



' is the contribution to the chemical 
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FIG. 1. Adsorbate structures calculated using DFT-GGA. 
(For the first eight diagrams analogous calculations were also 
performed for O in fee sites.) Small circles represent O atoms 
and large circles Ru atoms. 
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potential of the adsorbate due to the lateral interac- 
tions in the hamiltonian (Eq. 1) and is calculated here 
using transfer matrix techniques. Regarding the sticking 
coefficient, we note that dissociation is not activated ini- 
tially, but at (local) ctaierages of 6'> 0.5, it is hindered 
by an energy barrieralla. Under these circumstances, 
to obtain the coverage and temperature dependence of 
the sticking coefficient ab initio would be a significant 
undertaking. Therefore we use an analytic expression 
which well approximates the measured behavior in the 
temperature regime of desorptiontj: The sticking coef- 
ficient drops approximately as (2/3 — 6*) ^ and for cover- 
ages above 2/3 it remains very small up to a monolayer. 
The actual equation we use is: S'dis(^) = Soexp[—{9/(T)'^], 
with 5*0=0.27 and cr=0.3. 
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FIG. 2. Theoretical (left panel) and experimentalcj (right 
panel) TPD spectra (heating rate of 6 K/s). For the theoret- 
ical results initial coverages are O =0.1 to 1.0 in steps of 0.1; 
the experimental results also span the initial coverage region 
of 6* to IML. 

Our calculated TPD spectrail are shown in-Fig. 2 
which are compared to recent experimental dataEil. For 
low initial coverage we note that in the theory the oxygen 
desorbs at about 100 K higher temperature than in ex- 
periment, reflecting an overbinding of the O atoms. We 
believe that this size of error is typical for present day 
state-of-the-art calculations and don't see that (or how) 
better accuracy will be achieved in the coming five to ten 
years. Aside from this, it can be seen that the theoret- 
ical spectra exhibit all the features of the experimental 
data: (i) a shift of the peak maxima to lower tempera- 
tures for higher initial coverages (appropriate for second 
order desorption and/or for repulsive interactions) and 
(ii) a steepening of the leading edge for higher initial 
coverages. In our theory this steepening is a reflection of 
two facts: Firstly, for increasing coverages, the repulsive 



Site 
hep 
fee 



TABLE I. Adsorption energies (in eV) for O on Ru (0001) 
with respect to I/2O2 for various coverages. 



next nearest neighbor interactions lower the binding 
energy and thus lower the onset of desorption, broaden- 
ing the TPD spectra and steepening the rising edge; sec- 
ondly, and more importantly, is the effect of the rapidly 
decreasing sticking probability for increasing coverageEj. 
Because the sticking coefficient is much smaller for cover- 
ages above 2/3ML, desorption is delayed to higher tem- 
peratures and the last third of a monolayer desorbs over 
a very narrow temperature range. If we use the ideal 
dissociative sticking, S{9) — 6*0(1 — 0)^, a similar steep- 
ening still occurs but the onset of desorption occurs 100 K 
earlier at a mosplayer. A similar behavior has been dis- 
cussed in Ref.E^. We can trace the two shoulders (at 
HOOK and 1300K) to the synergy of the interactions that 
at lower temperatures lead to the formation of the (2 x 1) 
and (2 x 2) ordered structures which will also be seen in 
the heat of adaaxption. 

Piercy et ali3 find that to obtain a satisfactory ex- 
planation of the surface phase diagram in the vicinity 
of the order-disorder transition temperature, a spillover 
into fee sites of about 12 % takes place. In the tem- 
perature regime of desorption, localization in hep sites 
should be even less. We find, however, that the overall 
features of the TPD spectra remain essentially unchanged 
whether spillover is included or not. We have also tested 
the importance of the trio interactions on desorption. Ne- 
glecting them increases the overall repulsion for coverages 
larger than 2/3 ML and consequently broadens the TPD 
spectra, reducing the agreement with experiment. We 
therefore conclude that for high 0-coverages trio interac- 
tions play an important role. 
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FIG. 3. The heat of adsorption as a function of coverage 
for temperatures of T=700 K (sharpest features), 1100 K, and 
1500 K. 
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Having the chemical potential as a function of temper- 
ature and coverage, we can also calculate the equilibrium 
properties of the adsorbate, such as adsorption isobars. 
As an example, we present in Fig. 3 the isosteric heat 
of adsorptionS for a few temperatures. At the highest 
temperature a smooth decrease is observed. At the low- 
est temperature, sharp peaks (and decreases) at 1/4, 1/2, 
3/4, and 1 ML occur. These coverages correspond to each 
oipithc orderftd phases that i)rm in nature— f..e. (2 x 2)- 
(M, (2 X l)-Oa, (2 X 2)-30Ll"y, and (1 x l)-^. (Note that 
the existence of the (2 x 2)-30, and (1 x l)-0 adsorption, 
phases was at first predicted by DFT-GGA calculationsEll 
and subsequently confirmed by experiment; a nice suc- 
cess of DFT.) The rises in between originate from the 
third neighbor attractions, and also from the trios for 
the higher coverages. Our results show no tendency for 
the stability of a 61 = 1/3 phase, i.e. a (\/3 x -\/3)i?30° 
structure, in agreement with experiment. We now turn 
to the surface phase diagram. From Tab. II it can be seen 
that the overall agreement of the interaction parameters 
determined from our density- functional jcalculations and 
from the best-fit to experiment in Ref.tj is, in general, 
astonishingly good; but there are some significant differ- 
ences (i.e. more than 50% for Vj[;'^) which could be ex- 
pected. We find that our Monte Carlo simulations yields 
a surface phase diagram rather similar to that of Piercy 
et oLEj where coverage up to half a monolayer was con- 
sidered. Our present simulations also included higher O 
coverages and correctly predict formation of the (2 x 2)- 
30 phase. 

In summary, we have presented a first-principles based 
approach for calculation of the thermodynamics and ki- 
netics of an adsorbate on a surface. We used density- 
functional theory to create a lattice gas hamiltonian from 
which we evaluated the partition function. Our theoret- 
ical temperature programmed thermal desorption spec- 
tra, heats of adsorption, and the surface phase diagram 
for O on Ru (0001) show very good overall agreement 
with available experimental results, providing confidence 
in our approach. We found that trio interactions, as well 
as the sticking coefficient, play an important role in the 
TPD spectra of the present system. The attractive trio 
interactions also apparently help stabilize the higher cov- 
erage (2 X l)-0, (2 X 2)-30, and (1 x l)-0 phases. The 
affect of spillover into fee sites was found to have a min- 
imal affect on the TPD spectra. 
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